Figures for Seasonally variable thermal performance curves prevent adverse effects of heatwaves

# Organization notes
# Main text includes effect size figures 

# Supp material includes: 
# Raw data plots for all four traits 
# Effect size plot showing all traits (done)
# Contrasts from modeling approach 

Main Text Figures

comb_epr_plot = comb_preds %>% 
  filter(metric == "EPR") %>% 
  ggplot() +
  geom_point(aes(temp, rate, colour = curve_id), 
             filter(comb_d, metric == "EPR"), 
             size = 1.5, alpha = 0.6, 
             position = position_jitter(width = 0.5, height = 0)) +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper, group = curve_id), 
              filter(comb_boot_conf_preds, metric == "EPR"), 
              fill = 'grey60', alpha = 0.3) +
  geom_line(aes(temp, .fitted, col = curve_id), linewidth = 2) +
  scale_colour_manual(values = comb_colors) + 
  labs(x = "", 
       y = "Egg Production \n(eggs/female/day)",
       colour = "Month") + 
  theme_matt(base_size = 12)

comb_hs_plot = comb_preds %>% 
  filter(metric == "HF") %>% 
  ggplot() +
  geom_point(aes(temp, rate, colour = curve_id), filter(comb_d, metric == "HF"), size = 1.5, alpha = 0.6, 
             position = position_jitter(width = 0.5, height = 0)) +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper, group = curve_id), filter(comb_boot_conf_preds, metric == "HF"), fill = 'grey60', alpha = 0.3) +
  geom_line(aes(temp, .fitted, col = curve_id), linewidth = 2) +
  #scale_color_brewer(type = "div", palette = 5, direction = -1) + 
  #scale_color_viridis_d(option = "mako") + 
  scale_colour_manual(values = comb_colors) + 
  labs(x = "", 
       y = "Hatching Success \n(%)",
       colour = "Month") + 
  theme_matt(base_size = 12)

comb_rf_plot = comb_preds %>% 
  filter(metric == "RF") %>% 
  ggplot() +
  geom_point(aes(temp, rate, colour = curve_id), filter(comb_d, metric == "RF"), size = 1.5, alpha = 0.6, 
             position = position_jitter(width = 0.5, height = 0)) +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper, group = curve_id), filter(comb_boot_conf_preds, metric == "RF"), fill = 'grey60', alpha = 0.3) +
  geom_line(aes(temp, .fitted, col = curve_id), linewidth = 2) +
  #scale_color_brewer(type = "div", palette = 5, direction = -1) + 
  #scale_color_viridis_d(option = "mako") + 
  scale_colour_manual(values = comb_colors) + 
  labs(x = "Temperature (°C)", 
       y = "Offspring Production \n(offspring/female/day)",
       colour = "Month") + 
  theme_matt(base_size = 12)

comb_tsc = ggplot(comb_surv, aes(x=Temp, y=Surv, colour=Month)) + 
  geom_point(size=1.5, position=position_jitter(width=0.1, height=0.03)) +
  xlab("Temperature (°C)")+
  ylab("Survivorship \n(proportion survived)")+
  labs(colour = "Month") + 
  geom_hline(yintercept = 0.5, linetype = "dashed") +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se=T, linewidth = 2) +
  scale_y_continuous(breaks = c(0,1)) + 
  #scale_color_brewer(type = "div", palette = 5, direction = -1) + 
  #scale_color_viridis_d(option = "mako") + 
  scale_colour_manual(values = comb_colors) + 
  theme_matt(base_size = 12)

ggarrange(comb_epr_plot, comb_hs_plot, comb_rf_plot, comb_tsc, 
          ncol = 2, nrow = 2,
          common.legend = T, legend = "bottom", 
          labels = "AUTO", vjust = 1)

# ggarrange(comb_epr_plot, comb_hs_plot, comb_rf_plot, comb_tsc, nrow = 1,
#           common.legend = T, legend = "bottom")

combined_opt_coll = comb_params %>% 
  filter(metric == "RF" & term == "topt") %>% 
  ggplot(aes(x = growth_temp, y = estimate, shape = species)) + 
  geom_smooth(data = filter(comb_params, metric == "RF" & term == "topt" & curve_id != "November_2"),
              method = "lm", colour = "grey50") + 
  geom_point(size = 3, stroke = 1) + 
  scale_shape_manual(values = c(16,21)) + 
  ylab("Optimum (°C)") + 
  xlab("Collection Temperature (°C)") + 
  labs(colour = "Month") + 
  theme_matt(base_size = 12)

combined_opt_diff = comb_params %>% 
  filter(metric == "RF" & term == "topt") %>% 
  ggplot(aes(x = growth_temp, y = margin, shape = species)) +
  geom_hline(yintercept = 0, linewidth =1, linetype = "dashed") +
  geom_smooth(data = filter(comb_params, metric == "RF" & term == "topt" & curve_id != "November_2"),
              method = "lm", colour = "grey50") + 
  geom_point(size = 3, stroke = 1) + 
  ylab("Margin (°C)") + 
  xlab("Collection Temperature (°C)") + 
  scale_shape_manual(values = c(16,21)) + 
  theme_matt(base_size = 12) 

combined_ld_coll = ggplot(combined_tolerance, aes(x = Coll_temp, y = LD50, shape = species)) + 
  geom_smooth(method = "lm", colour = "grey50") + 
  geom_point(size = 3, stroke = 1) + 
  scale_shape_manual(values = c(16,21)) + 
  ylab("Thermal Tolerance (°C)") + 
  xlab("Collection Temperature (°C)") + 
  theme_matt(base_size = 12)

combined_ld_diff = ggplot(combined_tolerance, aes(x = Coll_temp, y = margin, shape = species)) +
  geom_smooth(method = "lm", colour = "grey50") + 
  geom_point(size = 3, stroke = 1) + 
  scale_shape_manual(values = c(16,21)) + 
  ylab("Warming Tolerance (°C)") + 
  xlab("Collection Temperature (°C)") + 
  theme_matt(base_size = 12)

ggarrange(combined_opt_coll, combined_opt_diff, combined_ld_coll, combined_ld_diff, ncol = 2, nrow = 2, common.legend = T,
          legend = "bottom", labels = "AUTO")

f0_model_females = F0_epr %>% 
  group_by(Month, Treatment, Female) %>% 
  count() %>% 
  filter(n == 2) %>% 
  mutate('female_ID' = paste(Month, Treatment, Female, sep = "_"))

f0_model_data = F0_epr %>% 
  mutate('female_ID' = paste(Month, Treatment, Female, sep = "_")) %>% 
  filter(female_ID %in% f0_model_females$female_ID) %>% 
  mutate(Day = if_else(Day == "1_to_3", "Short", "Long"),
         Month = fct_relevel(Month, "June", "August", "November"),
         Day = fct_relevel(Day, "Short", "Long"))

F0_hatched.model = lme4::glmer(data = f0_model_data, family = poisson, 
                               Hatched ~ Treatment * Day * Month + (1|female_ID))

summary(F0_hatched.model)
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
##  Family: poisson  ( log )
## Formula: Hatched ~ Treatment * Day * Month + (1 | female_ID)
##    Data: f0_model_data
## 
##      AIC      BIC   logLik deviance df.resid 
##   7322.9   7375.8  -3648.5   7296.9      419 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -9.044 -1.385 -0.256  1.178  8.188 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  female_ID (Intercept) 1.427    1.195   
## Number of obs: 432, groups:  female_ID, 216
## 
## Fixed effects:
##                                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                              4.06299    0.19477  20.861  < 2e-16 ***
## TreatmentHeatwave                        0.16319    0.27815   0.587  0.55742    
## DayLong                                 -0.20893    0.02343  -8.918  < 2e-16 ***
## MonthAugust                              0.71281    0.27516   2.591  0.00958 ** 
## MonthNovember                           -4.73525    0.35172 -13.463  < 2e-16 ***
## TreatmentHeatwave:DayLong               -0.40687    0.03551 -11.457  < 2e-16 ***
## TreatmentHeatwave:MonthAugust           -0.08077    0.39427  -0.205  0.83769    
## TreatmentHeatwave:MonthNovember          0.69060    0.46905   1.472  0.14093    
## DayLong:MonthAugust                     -0.40799    0.03391 -12.030  < 2e-16 ***
## DayLong:MonthNovember                    2.32910    0.20062  11.609  < 2e-16 ***
## TreatmentHeatwave:DayLong:MonthAugust    0.32491    0.04852   6.697 2.13e-11 ***
## TreatmentHeatwave:DayLong:MonthNovember -0.06312    0.23667  -0.267  0.78968    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnH DayLng MnthAg MnthNv TrH:DL TrH:MA TrH:MN DyL:MA DyL:MN TH:DL:MA
## TretmntHtwv -0.698                                                                        
## DayLong     -0.054  0.038                                                                 
## MonthAugust -0.708  0.494  0.038                                                          
## MonthNovmbr -0.551  0.386  0.030  0.390                                                   
## TrtmntHt:DL  0.036 -0.050 -0.660 -0.025 -0.020                                            
## TrtmntHt:MA  0.494 -0.706 -0.027 -0.698 -0.272  0.035                                     
## TrtmntHt:MN  0.414 -0.593 -0.022 -0.293 -0.750  0.030  0.418                              
## DyLng:MnthA  0.037 -0.026 -0.691 -0.049 -0.021  0.456  0.034  0.015                       
## DyLng:MnthN  0.006 -0.004 -0.117 -0.004 -0.506  0.077  0.003  0.379  0.081                
## TrtmH:DL:MA -0.026  0.037  0.483  0.034  0.014 -0.732 -0.045 -0.022 -0.699 -0.056         
## TrtmH:DL:MN -0.005  0.008  0.099  0.004  0.429 -0.150 -0.005 -0.438 -0.068 -0.848  0.110

knitr::kable(car::Anova(F0_hatched.model, type = "III"))

	Chisq	Df	Pr(>Chisq)
(Intercept)	435.1722039	1	0.0000000
Treatment	0.3441978	1	0.5574162
Day	79.5279360	1	0.0000000
Month	253.7336545	2	0.0000000
Treatment:Day	131.2551698	1	0.0000000
Treatment:Month	2.9847250	2	0.2248408
Day:Month	304.0138675	2	0.0000000
Treatment:Day:Month	45.8612180	2	0.0000000

knitr::kable(emmeans::emmeans(F0_hatched.model, ~ Day | Treatment * Month, type = "response") %>% pairs())

contrast	Treatment	Month	ratio	SE	df	null	z.ratio	p.value
Short / Long	Control	June	1.2323565	0.0288718	Inf	1	8.917844	0
Short / Long	Heatwave	June	1.8511273	0.0494058	Inf	1	23.072465	0
Short / Long	Control	August	1.8532162	0.0454425	Inf	1	25.159055	0
Short / Long	Heatwave	August	2.0115022	0.0445963	Inf	1	31.522871	0
Short / Long	Control	November	0.1200108	0.0239124	Inf	1	-10.640644	0
Short / Long	Heatwave	November	0.1920146	0.0235545	Inf	1	-13.452169	0

knitr::kable(emmeans::emmeans(F0_hatched.model, ~ Treatment | Day * Month, type = "response") %>% pairs())

contrast	Day	Month	ratio	SE	df	null	z.ratio	p.value
Control / Heatwave	Short	June	0.8494307	0.2362733	Inf	1	-0.5866837	0.5574162
Control / Heatwave	Long	June	1.2759329	0.3555254	Inf	1	0.8745263	0.3818317
Control / Heatwave	Short	August	0.9208831	0.2572647	Inf	1	-0.2950316	0.7679697
Control / Heatwave	Long	August	0.9995371	0.2798544	Inf	1	-0.0016538	0.9986804
Control / Heatwave	Short	November	0.4257991	0.1608150	Inf	1	-2.2606223	0.0237827
Control / Heatwave	Long	November	0.6812692	0.2167769	Inf	1	-1.2061690	0.2277523

f0_effects = emmeans::emmeans(F0_hatched.model, ~ Treatment | Day * Month, type = "response") %>% pairs() %>% data.frame() %>% 
  mutate(Day = fct_relevel(Day, "Short", "Long"),
         Month = fct_relevel(Month, "June", "August", "November")) %>% 
  ggplot(aes(x = Day, y = ratio, shape = Day)) + 
  facet_wrap(.~Month, nrow = 1) + 
  geom_hline(yintercept = 1) + 
  geom_point(position = position_dodge(width = 0.3), 
             size = 7) + 
  geom_errorbar(aes(ymin = ratio - SE, ymax = ratio + SE),
                position = position_dodge(width = 0.3),
                width = 0.2, linewidth = 1) + 
  labs(x = "Duration", 
       y = "Ratio (Control:Heatwave)",
       shape = "Duration") + 
  theme_bw(base_size = 20) + 
  theme(legend.position = "none", 
        panel.grid = element_blank())

f1_model_data = F1_epr %>% 
  mutate(Offspring_temp = as.factor(Offspring_temp)) %>% 
  ungroup() %>% 
  mutate(Day = if_else(Day == "1_to_3", "Short", "Long"),
         Month = fct_relevel(Month, "June", "August", "November"),
         Day = fct_relevel(Day, "Short", "Long"))

f0_epr_plot = ggplot(f0_model_data, 
                 aes(x = Day, y = Total, colour = Treatment)) + 
  facet_grid(.~Month, scales = "free_y") + 
  geom_hline(yintercept = 0) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  labs(y = "Egg Production \n(per female per day)") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        axis.text.x = element_text(angle = 300, hjust = 0, vjust = 0.5),
        axis.title.x = element_blank(),
        strip.text = element_text(size = 18))

f1_epr_plot = ggplot(f1_model_data, 
                 aes(x = factor(Offspring_temp), y = Total, colour = Parental_treatment)) + 
  facet_grid(Month~Day) + 
  #geom_violin(position = position_dodge(width = 1)) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  scale_y_continuous(breaks = c(0, 100, 200)) + 
  labs(x = "Offspring Temp. (°C)", 
       y = "Egg Production \n(per female per day)",
       colour = "Parental Treatment") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        strip.text = element_text(size = 18))

raw_epr = ggarrange(f0_epr_plot, f1_epr_plot, nrow = 2, 
          heights = c(0.3, 0.7), common.legend = T, legend = "bottom")

###

f0_hs_plot = ggplot(f0_model_data, 
                 aes(x = Day, y = Success, colour = Treatment)) + 
  facet_grid(.~Month, scales = "free_y") + 
  geom_hline(yintercept = 0) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  labs(y = "Hatching Success") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        axis.text.x = element_text(angle = 300, hjust = 0, vjust = 0.5),
        axis.title.x = element_blank(),
        strip.text = element_text(size = 18))

f1_hs_plot = ggplot(f1_model_data, 
                 aes(x = factor(Offspring_temp), y = Success, colour = Parental_treatment)) + 
  facet_grid(Month~Day) + 
  #geom_violin(position = position_dodge(width = 1)) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  scale_y_continuous(breaks = c(0, 1)) + 
  labs(x = "Offspring Temp. (°C)", 
       y = "Hatching Success",
       colour = "Parental Treatment") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        strip.text = element_text(size = 18))

raw_hs = ggarrange(f0_hs_plot, f1_hs_plot, nrow = 2, 
          heights = c(0.3, 0.7), common.legend = T, legend = "bottom")

###

f0_prod_plot = ggplot(f0_model_data, 
                 aes(x = Day, y = Hatched, colour = Treatment)) + 
  facet_grid(.~Month, scales = "free_y") + 
  geom_hline(yintercept = 0) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  labs(y = "Offspring Production \n(per female per day)") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        axis.text.x = element_text(angle = 300, hjust = 0, vjust = 0.5),
        axis.title.x = element_blank(),
        strip.text = element_text(size = 18))

f1_prod_plot = ggplot(f1_model_data, 
                 aes(x = factor(Offspring_temp), y = Hatched, colour = Parental_treatment)) + 
  facet_grid(Month~Day) + 
  #geom_violin(position = position_dodge(width = 1)) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  scale_y_continuous(breaks = c(0, 100, 200)) + 
  labs(x = "Offspring Temp. (°C)", 
       y = "Offspring Production \n(per female per day)",
       colour = "Parental Treatment") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        strip.text = element_text(size = 18))

raw_production = ggarrange(f0_prod_plot, f1_prod_plot, nrow = 2, 
          heights = c(0.3, 0.7), common.legend = T, legend = "bottom")

ggarrange(raw_epr, raw_hs, raw_production, nrow = 1, common.legend = T, legend = "bottom",
          labels = "AUTO")

# F1_hatched.model = lme4::glmer(data = f1_model_data, family = poisson, 
#                                Hatched ~ Parental_treatment * Day * Offspring_temp + 
#                                  (1 + Parental_treatment|Month))

F1_hatched.model = glm(data = f1_model_data, family = "poisson",
                       Hatched ~ Parental_treatment * Day * Offspring_temp * Month)

summary(F1_hatched.model)
## 
## Call:
## glm(formula = Hatched ~ Parental_treatment * Day * Offspring_temp * 
##     Month, family = "poisson", data = f1_model_data)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -15.6120   -4.1732   -0.6726    1.8966   22.9798  
## 
## Coefficients: (1 not defined because of singularities)
##                                                                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                                        3.312e+00  3.486e-02  95.008  < 2e-16 ***
## Parental_treatmentHeatwave                                        -3.379e-01  5.402e-02  -6.255 3.97e-10 ***
## DayLong                                                           -1.554e+00  8.344e-02 -18.623  < 2e-16 ***
## Offspring_temp17                                                   1.181e+00  3.985e-02  29.635  < 2e-16 ***
## Offspring_temp22                                                   1.209e+00  3.988e-02  30.304  < 2e-16 ***
## MonthAugust                                                       -1.122e+00  7.033e-02 -15.953  < 2e-16 ***
## MonthNovember                                                     -1.968e+00  1.515e-01 -12.990  < 2e-16 ***
## Parental_treatmentHeatwave:DayLong                                -6.622e-01  1.559e-01  -4.249 2.15e-05 ***
## Parental_treatmentHeatwave:Offspring_temp17                        8.120e-01  5.936e-02  13.679  < 2e-16 ***
## Parental_treatmentHeatwave:Offspring_temp22                        3.063e-01  6.068e-02   5.048 4.47e-07 ***
## DayLong:Offspring_temp17                                           1.864e+00  8.723e-02  21.370  < 2e-16 ***
## DayLong:Offspring_temp22                                           1.517e+00  8.798e-02  17.237  < 2e-16 ***
## Parental_treatmentHeatwave:MonthAugust                            -1.351e-01  1.124e-01  -1.201 0.229590    
## Parental_treatmentHeatwave:MonthNovember                           1.635e+01  8.570e+01   0.191 0.848711    
## DayLong:MonthAugust                                                2.054e+00  1.138e-01  18.047  < 2e-16 ***
## DayLong:MonthNovember                                             -2.498e+00  1.014e+00  -2.463 0.013786 *  
## Offspring_temp17:MonthAugust                                      -4.293e-01  9.084e-02  -4.726 2.30e-06 ***
## Offspring_temp22:MonthAugust                                      -6.995e-01  8.742e-02  -8.002 1.23e-15 ***
## Offspring_temp17:MonthNovember                                    -3.275e-01  1.856e-01  -1.765 0.077602 .  
## Offspring_temp22:MonthNovember                                    -1.585e+01  8.569e+01  -0.185 0.853203    
## Parental_treatmentHeatwave:DayLong:Offspring_temp17               -5.423e-01  1.604e-01  -3.380 0.000724 ***
## Parental_treatmentHeatwave:DayLong:Offspring_temp22                3.778e-01  1.612e-01   2.344 0.019102 *  
## Parental_treatmentHeatwave:DayLong:MonthAugust                    -2.488e-04  2.081e-01  -0.001 0.999046    
## Parental_treatmentHeatwave:DayLong:MonthNovember                  -1.434e+01  8.569e+01  -0.167 0.867053    
## Parental_treatmentHeatwave:Offspring_temp17:MonthAugust           -1.349e-01  1.327e-01  -1.017 0.309325    
## Parental_treatmentHeatwave:Offspring_temp22:MonthAugust            5.513e-01  1.313e-01   4.198 2.69e-05 ***
## Parental_treatmentHeatwave:Offspring_temp17:MonthNovember         -1.503e+01  8.570e+01  -0.175 0.860799    
## Parental_treatmentHeatwave:Offspring_temp22:MonthNovember          1.006e+00  1.369e+00   0.735 0.462321    
## DayLong:Offspring_temp17:MonthAugust                              -7.350e-01  1.299e-01  -5.657 1.54e-08 ***
## DayLong:Offspring_temp22:MonthAugust                              -9.589e-01  1.301e-01  -7.372 1.69e-13 ***
## DayLong:Offspring_temp17:MonthNovember                            -2.574e+00  1.241e+00  -2.074 0.038055 *  
## DayLong:Offspring_temp22:MonthNovember                             1.380e+01  8.569e+01   0.161 0.872056    
## Parental_treatmentHeatwave:DayLong:Offspring_temp17:MonthAugust    6.825e-02  2.244e-01   0.304 0.761005    
## Parental_treatmentHeatwave:DayLong:Offspring_temp22:MonthAugust   -7.590e-01  2.280e-01  -3.329 0.000872 ***
## Parental_treatmentHeatwave:DayLong:Offspring_temp17:MonthNovember  2.040e+01  8.569e+01   0.238 0.811813    
## Parental_treatmentHeatwave:DayLong:Offspring_temp22:MonthNovember         NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 62808  on 930  degrees of freedom
## Residual deviance: 25570  on 896  degrees of freedom
## AIC: 28981
## 
## Number of Fisher Scoring iterations: 11
#car::Anova(F1_hatched.model,test = "F", type = "III")

f1_hatched_contrasts = emmeans::emmeans(F1_hatched.model, ~ Parental_treatment | Day * Month * Offspring_temp, type = "response") %>% pairs() %>% data.frame() %>% 
  #mutate(estimate = estimate * -1) %>% 
  drop_na() %>% 
  mutate("ID" = paste0(Month, Offspring_temp))

f1_effects = ggplot(f1_hatched_contrasts, aes(x = Day, y = ratio, group = ID)) + 
  facet_grid(Offspring_temp~Month) + 
  geom_hline(yintercept = 1) +
  geom_line(linewidth = 0.8) + 
  geom_point(size = 2) + 
  geom_errorbar(aes(ymin = ratio - SE, ymax = ratio + SE), width = 0.1, linewidth = 1) + 
  labs(x = "Duration", 
       y = "Ratio (Control:Heatwave)") + 
  theme_bw(base_size = 20) + 
  theme(panel.grid = element_blank())

ggarrange(f0_effects, f1_effects, nrow = 2, heights = c(0.3, 0.7))

f1_size_data = F1_fbs %>% 
  ungroup() %>% 
  mutate(Day = if_else(Day == "1_to_3", "Short", "Long"),
         Month = fct_relevel(Month, "June", "August", "November"),
         Day = fct_relevel(Day, "Short", "Long"))

F1_size.model = glm(data = f1_size_data,
                    Size ~ Parental_treatment * Day * Offspring_temp * Month)

summary(F1_size.model)
## 
## Call:
## glm(formula = Size ~ Parental_treatment * Day * Offspring_temp * 
##     Month, data = f1_size_data)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.145900  -0.025187  -0.002125   0.027106   0.156603  
## 
## Coefficients:
##                                                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                                                      1.1129676  0.0191453  58.133  < 2e-16 ***
## Parental_treatmentHeatwave                                      -0.2021937  0.0314027  -6.439 2.02e-10 ***
## DayLong                                                         -0.0639884  0.0270852  -2.362 0.018380 *  
## Offspring_temp                                                  -0.0116568  0.0010991 -10.606  < 2e-16 ***
## MonthAugust                                                     -0.0778029  0.0286411  -2.716 0.006734 ** 
## MonthNovember                                                   -0.1583042  0.0549167  -2.883 0.004044 ** 
## Parental_treatmentHeatwave:DayLong                               0.1530333  0.0414772   3.690 0.000239 ***
## Parental_treatmentHeatwave:Offspring_temp                        0.0105033  0.0017407   6.034 2.40e-09 ***
## DayLong:Offspring_temp                                           0.0033498  0.0015625   2.144 0.032333 *  
## Parental_treatmentHeatwave:MonthAugust                           0.1597975  0.0426050   3.751 0.000188 ***
## Parental_treatmentHeatwave:MonthNovember                         0.0138210  0.0735478   0.188 0.850986    
## DayLong:MonthAugust                                             -0.0208334  0.0394776  -0.528 0.597827    
## DayLong:MonthNovember                                            0.1684485  0.0627450   2.685 0.007404 ** 
## Offspring_temp:MonthAugust                                       0.0030094  0.0016350   1.841 0.066035 .  
## Offspring_temp:MonthNovember                                     0.0104002  0.0037188   2.797 0.005281 ** 
## Parental_treatmentHeatwave:DayLong:Offspring_temp               -0.0074358  0.0023391  -3.179 0.001533 ** 
## Parental_treatmentHeatwave:DayLong:MonthAugust                  -0.0655264  0.0576266  -1.137 0.255826    
## Parental_treatmentHeatwave:DayLong:MonthNovember                 0.0366208  0.0858275   0.427 0.669723    
## Parental_treatmentHeatwave:Offspring_temp:MonthAugust           -0.0104701  0.0023926  -4.376 1.36e-05 ***
## Parental_treatmentHeatwave:Offspring_temp:MonthNovember         -0.0008553  0.0044990  -0.190 0.849261    
## DayLong:Offspring_temp:MonthAugust                               0.0021581  0.0022665   0.952 0.341278    
## DayLong:Offspring_temp:MonthNovember                            -0.0099254  0.0040868  -2.429 0.015364 *  
## Parental_treatmentHeatwave:DayLong:Offspring_temp:MonthAugust    0.0029113  0.0032879   0.885 0.376166    
## Parental_treatmentHeatwave:DayLong:Offspring_temp:MonthNovember -0.0041042  0.0051223  -0.801 0.423219    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.001721073)
## 
##     Null deviance: 2.7198  on 863  degrees of freedom
## Residual deviance: 1.4457  on 840  degrees of freedom
##   (2 observations deleted due to missingness)
## AIC: -3021.6
## 
## Number of Fisher Scoring iterations: 2
car::Anova(F1_size.model, type = "III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: Size
##                                             LR Chisq Df Pr(>Chisq)    
## Parental_treatment                            41.457  1  1.205e-10 ***
## Day                                            5.581  1  0.0181527 *  
## Offspring_temp                               112.484  1  < 2.2e-16 ***
## Month                                         12.730  2  0.0017203 ** 
## Parental_treatment:Day                        13.613  1  0.0002246 ***
## Parental_treatment:Offspring_temp             36.410  1  1.599e-09 ***
## Day:Offspring_temp                             4.596  1  0.0320460 *  
## Parental_treatment:Month                      15.161  2  0.0005104 ***
## Day:Month                                      9.125  2  0.0104334 *  
## Offspring_temp:Month                           9.540  2  0.0084791 ** 
## Parental_treatment:Day:Offspring_temp         10.105  1  0.0014784 ** 
## Parental_treatment:Day:Month                   2.062  2  0.3566473    
## Parental_treatment:Offspring_temp:Month       20.328  2  3.852e-05 ***
## Day:Offspring_temp:Month                       8.623  2  0.0134123 *  
## Parental_treatment:Day:Offspring_temp:Month    2.110  2  0.3482539    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

day_cols = c("Short" = "grey80", "Long" = "grey30")

size_temp1 = ggplot(f1_size_data, aes(x = Offspring_temp, y = Size, colour = Parental_treatment)) + 
  facet_grid(Month~Day) + 
  geom_jitter(width = 0.5, size = 1.6, alpha = 0.4) +
  geom_smooth(method = "lm", linewidth = 1.4, alpha = 0.2) + 
  labs(x = "Offspring Temperature (°C)",
       y = "Size (mm)",
       colour = "Parental Treatment") + 
  scale_x_continuous(breaks = c(12,17,22)) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  theme_matt(base_size = 15) + theme(legend.position = "bottom",
                                     panel.grid = element_blank(),
                                     panel.border = element_rect(fill = NA, colour = "black"))

#car::Anova(f1_size.model, type = "III")
size_temp2 = emmeans::emtrends(F1_size.model, c("Month", "Day", "Parental_treatment"), var = "Offspring_temp") %>% 
  as.data.frame() %>% 
  ggplot(aes(x = Parental_treatment, y = Offspring_temp.trend, 
             colour = Day, group = Day)) + 
  facet_wrap(Month~.) + 
  geom_hline(yintercept = 0) + 
  geom_line(linewidth = 1.5,
            position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                linewidth = 1.5, width = 0.5,
                position = position_dodge(width = 0.5)) + 
  geom_point(size = 3,
             position = position_dodge(width = 0.5)) + 
  scale_colour_manual(values = day_cols) + 
  labs(x = "Parental Treatment",  
       y = "Size Slope (mm / °C)") + 
  theme_matt(base_size = 15) + 
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 290,
                                   hjust = 0, vjust = 0.5),
        plot.background = element_rect(fill = "white"))

ggarrange(size_temp1, size_temp2, labels = "AUTO", common.legend = F, legend = "bottom")

f1_size_contrasts = emmeans::emmeans(F1_size.model, ~ Parental_treatment | Day * Month * Offspring_temp, type = "response", at = list(Offspring_temp = c(12,17,22))) %>% pairs() %>% data.frame() %>% 
  drop_na() %>% 
  mutate("ID" = paste0(Month, Offspring_temp)) %>% 
  filter(!(ID == "November12" & Day == "Short")) %>% 
  filter(!(ID == "November22" & Day == "Short"))

ggplot(f1_size_contrasts, aes(x = Day, y = estimate, group = ID)) + 
  facet_grid(Offspring_temp~Month) + 
  geom_hline(yintercept = 0) +
  geom_line(linewidth = 0.8) + 
  geom_point(size = 2) + 
  geom_errorbar(aes(ymin = estimate - SE, ymax = estimate + SE), width = 0.1, linewidth = 1) + 
  labs(x = "Duration", 
       y = "Effect (mm; Control - Heatwave)") + 
  theme_bw(base_size = 20) + 
  theme(panel.grid = element_blank())

# F0_rf_summary$month = factor(F0_rf_summary$month, levels = c("June", "August", "November"))
# F0_rf_summary$duration = factor(F0_rf_summary$duration, levels = c("short", "long"))
# F0_dur_summary$month = factor(F0_dur_summary$month, levels = c("June", "August", "November"))
# 
# param_list = list("colour" = "black",
#                   "width" = 0.2) 
# 
# RF_short_db = plot(F0_RF_short, 
#                    axes.title.fontsize = 10,
#                    tick.fontsize = 10,
#                    effsize.markersize = 3,
#                    swarmplot.params = param_list,
#                    rawplot.ylabel = "Production",
#                    theme = ggpubr::theme_pubr())
# 
# RF_long_db = plot(F0_RF_long, 
#                   axes.title.fontsize = 10,
#                   tick.fontsize = 10,
#                   effsize.markersize = 3,
#                   swarmplot.params = param_list,
#                   rawplot.ylabel = "Production",
#                   theme = ggpubr::theme_pubr())
# 
# b1 = ggplot() + theme_pubclean() + ggtitle("          Short Heat Waves")
# b2 = ggplot() + theme_pubclean() + ggtitle("          Long Heat Waves")
# F0_fecundity_plot = ggarrange(b1, RF_short_db, b2, RF_long_db, ncol = 1, nrow = 4, heights = c(0.1, 1, 0.1, 1))

# x.axis_labels = c("1" = "short", "2" = "long", "3" = "short", "4" = "long")
# 
# F0_grid = F0_rf_summary %>% 
#   mutate(month = fct_relevel(month, "June", "August", "November")) %>% 
#   ggplot(aes(x = duration, y = difference, colour = trait, shape = duration)) + 
#   facet_grid(. ~ month) + 
#   geom_hline(yintercept = 0, colour = "black", linewidth = 0.3) +
#   geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high), width = 0, linewidth = 1.3) + 
#   geom_point(size = 5, fill = "white") + 
#   scale_colour_manual(values = c("body size" = "grey75", "production" = "black")) + 
#   scale_shape_manual(values = c("long" = 16, "short" = 21)) + 
#   scale_x_discrete(labels= x.axis_labels) +
#   ggtitle("Direct Effects (F0)") + 
#   xlab("") +
#   ylab("Effect Size\nHeatwave - Control") + 
#   ylim(-1,1.1) + 
#   theme_bw(base_size = 12) + 
#   theme(panel.grid = element_blank(),
#         axis.text = element_text(colour = "black"),
#         axis.text.x = element_text(angle = -45, hjust = 0, vjust = 0.5),
#         legend.position = "none")
# 
# F1_summary = bind_rows(F1_rf_effect_size, F1_bs_effect_size) %>% 
#   dplyr::select(variable, difference, 
#                 bca_ci_low, bca_ci_high, 
#                 month, duration, trait, generation, off_temp) %>% 
#   mutate("order_code" = paste(trait, duration, sep = "_"),
#          "order_number" = case_when(
#            order_code == "production_short" ~ 1,
#            order_code == "production_long" ~ 2,
#            order_code == "body size_short" ~ 3,
#            order_code == "body size_long" ~ 4),
#          month = fct_relevel(month, "June", "August", "November"))
# 
# F1_summary$order_number = factor(F1_summary$order_number, levels = c(1,2,3,4))
# F1_grid = ggplot(F1_summary, aes(x = order_number, y = difference, colour = trait, shape = duration, group = trait)) +
#   facet_grid(off_temp ~ month, ) + 
#   geom_hline(yintercept = 0, colour = "black", linewidth = 0.3) +
#   geom_line() + 
#   geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high), width = 0, linewidth = 1) + 
#   geom_point(size = 3, fill = "white") + 
#   scale_colour_manual(values = c("body size" = "grey75", "production" = "black")) + 
#   scale_shape_manual(values = c("long" = 16, "short" = 21)) + 
#   xlim(0.5,4.5) + 
#   scale_x_discrete(labels= x.axis_labels) +
#   xlab("") +
#   ylab("Effect Size\nHeatwave - Control") + 
#   ggtitle("Transgenerational Effects (F1)") + 
#   theme_bw(base_size = 12) + 
#   theme(panel.grid = element_blank(),
#         strip.background.x = element_blank(),
#         strip.text.x = element_blank(),
#         axis.text = element_text(colour = "black"),
#         axis.text.x = element_text(angle = -45, hjust = 0, vjust = 0.5))
# 
# ggarrange(F0_grid, F1_grid, nrow = 2, ncol = 1, heights = c(0.45,1), common.legend = T, legend = "right")

Supplemental Information

#field tpc parameters
comb_params %>%  
  mutate(curve_id = fct_relevel(curve_id, c("January", "February", "March", "April", "May", "June", 
                                            "July", "August", "September", "October", "November_1", "November_2"))) %>% 
  ggplot(aes(x = curve_id, y = estimate, colour = species)) +
  facet_wrap(term~metric, scales = 'free_y') + 
  geom_point(size = 4) +
  geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
  scale_colour_manual(values = c("royalblue1", "indianred2")) + 
  labs(x = "Month",
       y = "Parameter Estimate",
       colour = "Species") + 
  theme_bw(base_size = 16) +
  theme(panel.grid = element_blank(),
        axis.text.x = element_text(angle = 315, hjust = 0, vjust = 0.5))

#field tpc parameters
ggplot(comb_params, aes(x = growth_temp, y = estimate, colour = species)) +
  facet_wrap(term~metric, scales = 'free_y') + 
  geom_smooth(data = filter(comb_params, curve_id != "November_2"),
              method = "lm", se = F) + 
  geom_point(size = 4) +
  geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
  scale_colour_manual(values = c("royalblue1", "indianred2")) + 
  labs(x = "Collection Temperature (°C)",
       y = "Parameter Estimate",
       colour = "Species") + 
  theme_bw(base_size = 12) +
  theme(panel.grid = element_blank())

comb_params %>% 
  filter(term == "topt") %>% 
  ggplot(aes(x = growth_temp, y = margin, colour = species, group = species)) + 
  facet_grid(metric~., scales = 'free_y') + 
  geom_hline(yintercept = 0) +
  geom_smooth(data = filter(comb_params, term == "topt" & curve_id != "November_2"),
              method = "lm", se = F, linewidth = 1, colour = "black") + 
  geom_point(size = 4) +
  scale_colour_manual(values = c("royalblue1", "indianred2")) + 
  labs(x = "Collection Temperature",
       y = "Safety Margin") + 
  theme_bw(base_size = 18) +
  theme(panel.grid = element_blank())

# #Subsequent Rows
# F1_hs_effect_size$trait = "success"
# F1_hs_effect_size$generation = "F1"
# 
# F1_total_effect_size$trait = "epr"
# F1_total_effect_size$generation = "F1"
# 
# F1_rf_effect_size$trait = "production"
# F1_rf_effect_size$generation = "F1"
# 
# F1_bs_effect_size$trait = "body size"
# F1_bs_effect_size$generation = "F1"
# 
# F0_data = bind_rows(F0_hs_summary, F0_total_summary,F0_rf_summary, F0_size_summary) %>% 
#   dplyr::select(trait, difference, bca_ci_low, bca_ci_high, month, duration, trait, generation) %>% 
#   mutate("order_code" = paste(trait, duration, sep = "_"),
#          "order_number" = case_when(
#            order_code == "total_short" ~ 1,
#            order_code == "total_long" ~ 2,
#            order_code == "success_short" ~ 3,
#            order_code == "success_long" ~ 4,
#            order_code == "production_short" ~ 5,
#            order_code == "production_long" ~ 6,
#            order_code == "size_long" ~ 7),
#          month = fct_relevel(month, "June", "August", "November"),
#          trait = fct_relevel(trait, "total", "success", "production", "size"),
#          duration = fct_relevel(duration, "short", "long"),
#          group_ID = paste(month, trait, sep = "_"))
# 
# F0_data$order_number = factor(F0_data$order_number, levels = c(1,2,3,4,5,6,7))
# 
# 
# F1_data = bind_rows(F1_total_effect_size, F1_hs_effect_size, F1_rf_effect_size, F1_bs_effect_size) %>% 
#   dplyr::select(trait, difference, bca_ci_low, bca_ci_high, month, duration, generation, off_temp) %>% 
#   mutate("order_code" = paste(trait, duration, sep = "_"),
#          "order_number" = case_when(
#            order_code == "epr_short" ~ 1,
#            order_code == "epr_long" ~ 2,
#            order_code == "success_short" ~ 3,
#            order_code == "success_long" ~ 4,
#            order_code == "production_short" ~ 5,
#            order_code == "production_long" ~ 6,
#            order_code == "body size_short" ~ 7,
#            order_code == "body size_long" ~ 8),
#          trait = if_else(trait == "epr", "total", trait),
#          month = fct_relevel(month, "June", "August", "November"),
#          trait = fct_relevel(trait, "total", "success", "production", "size"),
#          duration = fct_relevel(duration, "short", "long"))
# 
# 
# F1_data$order_number = factor(F1_data$order_number, levels = c(1,2,3,4,5,6,7,8))
# 
# #Top row - F0 (direct effects)
# x.axis_labels = c("1" = "short", "2" = "long", "3" = "short", "4" = "long", 
#                   "5" = "short", "6" = "long", "7" = "short", "8" = "long")
# 
# F0_grid = ggplot(F0_data, aes(x = duration, y = difference, colour = trait, group = group_ID)) + 
#   facet_grid(. ~ month) + 
#   geom_line(position = position_dodge(width = 0.7),
#             linewidth = 1) + 
#   geom_hline(yintercept = 0, colour = "black", linewidth = 0.3) +
#   geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high), width = 0, linewidth = 1, 
#                 position = position_dodge(width = 0.7)) + 
#   geom_point(size = 4, fill = "white", position = position_dodge(width = 0.7)) + 
#   scale_shape_manual(values = c("long" = 16, "short" = 21)) + 
#   scale_x_discrete(labels= x.axis_labels) +
#   scale_colour_manual(values = c("size" = "darkgrey",
#                                  "success" = "gold",
#                                  "production" = "forestgreen",
#                                  "total" = "cornflowerblue")) +  
#   xlab("") +
#   ylab("Effect Size\nHeatwave - Control") + 
#   theme_bw(base_size = 12) + 
#   theme(panel.grid = element_blank(),
#         axis.text = element_text(colour = "black"),
#         axis.text.x = element_text(angle = -45, hjust = 0, vjust = 0.5),
#         legend.position = "none")
# 
# #Following three rows - F1 (transgeneration / indirect effects)
# F1_grid = ggplot(F1_data, aes(x = duration, y = difference, colour = trait, group = trait)) + 
#   facet_grid(off_temp ~ month, ) + 
#   geom_hline(yintercept = 0, colour = "black", linewidth = 0.3) +
#   geom_line(position = position_dodge(width = 0.5),
#             linewidth = 1) + 
#   geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high), width = 0, linewidth = 1,
#                 position = position_dodge(width = 0.5)) + 
#   geom_point(size = 3, fill = "white", position = position_dodge(width = 0.5)) + 
#   scale_shape_manual(values = c("long" = 16, "short" = 21)) + 
#   xlim(0.5,4.5) + 
#   scale_x_discrete(labels= x.axis_labels) +
#   scale_colour_manual(values = c("size" = "darkgrey",
#                                  "success" = "gold",
#                                  "production" = "forestgreen",
#                                  "total" = "cornflowerblue")) + 
#   xlab("") +
#   ylab("Effect Size \nHeatwave - Control") + 
#   theme_bw(base_size = 12) + 
#   theme(panel.grid = element_blank(),
#         strip.background.x = element_blank(),
#         strip.text.x = element_blank(),
#         axis.text = element_text(colour = "black"),
#         axis.text.x = element_text(angle = -45, hjust = 0, vjust = 0.5))
# 
# ggarrange(F0_grid, F1_grid, nrow = 2, ncol = 1, heights = c(0.35,1), common.legend = T, legend = "right", labels = "AUTO")

# effect_corr = F1_summary %>% 
#   select(trait, difference, month, duration, off_temp) %>%  
#   pivot_wider(id_cols = c(month, duration, off_temp),
#               names_from = trait, 
#               values_from = difference)
# 
# ggplot(effect_corr, aes(x = `body size`, y = production)) + 
#   geom_hline(yintercept = 0) + 
#   geom_vline(xintercept = 0) + 
#   geom_point(size = 3) + 
#   geom_smooth(method = "lm", se = F,
#               colour = "grey60",
#               size = 1) + 
#   labs(x = "Body Size Effect",
#        y = "Production Effect") + 
#   theme_matt()

#Effect of heatwave duration WITHIN treatment

f0_model_data %>%  
  group_by(Treatment, Day, Month) %>% 
  summarise(mean_hatched = mean(Hatched, na.rm = T)) %>% 
  ungroup() %>% 
  pivot_wider(id_cols = c("Treatment", "Month"),
              values_from = mean_hatched, 
              names_from = Day) %>% 
  mutate("effect" = Long - Short)
## # A tibble: 6 × 5
##   Treatment Month      Short  Long effect
##   <chr>     <fct>      <dbl> <dbl>  <dbl>
## 1 Control   June     104.    84.4  -19.6 
## 2 Control   August   125.    67.2  -57.3 
## 3 Control   November   0.871  7.26   6.39
## 4 Heatwave  June     108.    58.3  -49.6 
## 5 Heatwave  August   170.    84.4  -85.4 
## 6 Heatwave  November   2.2   11.5    9.26

duration.model = lme4::lmer(data = f0_model_data, 
                            Hatched ~ Treatment * Day * Month + (1|Month))

duration_pairs = emmeans::emmeans(duration.model, 
                                  ~ Day | Month * Treatment) %>% 
  pairs()

as.data.frame(summary(duration_pairs))[c('Month', 
                                         'Treatment', 
                                         'contrast', 
                                         'estimate', 
                                         'SE')] %>% 
  mutate(Month = fct_relevel(Month, c("June", "August", "November")),
         estimate = estimate * -1) %>% #Flips sign to make contrast Long - Short
  ggplot(aes(x = Month, fill = Treatment, y = estimate)) + 
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), colour = "black", size = 1) + 
  geom_errorbar(aes(ymin = estimate - SE, ymax = estimate + SE),
                width = 0.1, linewidth = 1,
                position = position_dodge(width = 0.9)) + 
  geom_hline(yintercept = 0) + 
  scale_fill_manual(values = c("grey30", "white")) + 
  labs(x = "",
       y = "Duration Contrast \n Long - Short events") + 
  theme_pubr(base_size = 18)

# How does heat wave duration affect transgenerational effects? Reaction norms shown below for effect size comparisons (heatwave - control) for different duration of parental exposure
# F1_data %>% 
#   dplyr::select(month, duration, off_temp, trait, difference, bca_ci_low, bca_ci_high) %>% 
#   mutate("ID" = paste(month, off_temp, trait, sep = "_"),
#          month = fct_relevel(month, c("June", "August", "November")),
#          duration = fct_relevel(duration, c("short", "long"))) %>% 
#   ggplot(aes(x = duration, y = difference, colour = month, group = ID)) + 
#   facet_wrap(trait~.) + 
#   geom_hline(yintercept = 0) + 
#   geom_line(size = 1, position = position_dodge(width = 0.1)) + 
#   geom_point(size = 2, position = position_dodge(width = 0.1)) + 
#   geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high),
#                 size = 0.75, width = 0.1,
#                 position = position_dodge(width = 0.1)) + 
#   labs(x = "Parental Exposure Duration", 
#        y = "Effect Size (Hedge's g) \n Heatwave - Control") + 
#   ylim(-5,5) + 
#   theme_pubr(base_size = 18)

# #Pulls out reaction norms where there's a sign change (changes from positive, neutral, or negative between duration groups) 
# 
# duration_effects = F1_data %>% 
#   dplyr::select(month, duration, off_temp, trait, difference, bca_ci_low, bca_ci_high) %>% 
#   mutate("ID" = paste(month, off_temp, trait, sep = "_"),
#          month = fct_relevel(month, c("June", "August", "November")),
#          duration = fct_relevel(duration, c("short", "long"))) %>% 
#   group_by(ID, trait, duration) %>%  
#   mutate("zero_diff" = case_when(
#     sign(bca_ci_low) == sign(bca_ci_high) ~ "does not overlap zero",
#     sign(bca_ci_low) != sign(bca_ci_high) ~ "overlaps zero"
#   )) %>% 
#   ungroup(duration) %>% 
#   mutate("change" = case_when(
#     sign(difference)[1] == sign(difference)[2] & zero_diff[1] == zero_diff[2] ~ "Same",
#     sign(difference)[1] != sign(difference)[2] & zero_diff[1] == "overlaps zero" & zero_diff[2] == "overlaps zero" ~ "Same",
#     sign(difference)[1] != sign(difference)[2] & zero_diff[1] == zero_diff[2] & zero_diff[1] == "does not overlap zero" ~ "Different",
#     sign(difference)[1] != sign(difference)[2] & zero_diff[1] != zero_diff[2] ~ "Different",
#     sign(difference)[1] == sign(difference)[2] & zero_diff[1] != zero_diff[2] ~ "Different"
#   )) %>% 
#   arrange(month, off_temp, trait) %>% 
#   filter(change == "Different")
# 
# select_rnorms = duration_effects %>% 
#   dplyr::select(-duration, -difference, -bca_ci_low, -bca_ci_high, -ID, -zero_diff) %>% 
#   distinct()
# 
# sig_changes = F1_data %>% 
#   dplyr::select(month, duration, off_temp, trait, difference, bca_ci_low, bca_ci_high) %>% 
#   mutate("ID" = paste(month, off_temp, trait, sep = "_"),
#          month = fct_relevel(month, c("June", "August", "November")),
#          duration = fct_relevel(duration, c("short", "long"))) %>% 
#   filter(ID %in% select_rnorms$ID)
# 
# ggplot(sig_changes, aes(x = duration, y = difference, colour = month, group = ID)) + 
#   facet_wrap(trait~.) + 
#   geom_hline(yintercept = 0) + 
#   geom_line(size = 1, position = position_dodge(width = 0.1)) + 
#   geom_point(size = 2, position = position_dodge(width = 0.1)) + 
#   geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high),
#                 size = 0.75, width = 0.1,
#                 position = position_dodge(width = 0.1)) + 
#   geom_label_repel(data = sig_changes %>% filter(duration == "long"), 
#                    aes(label = off_temp, 
#                        x = duration,
#                        y = difference, 
#                        color = month),
#                    box.padding = 0.5,
#                    nudge_x = 0.2,
#                    size = 7,
#                    show.legend=FALSE) + 
#   labs(x = "Parental Exposure Duration", 
#        y = "Effect Size (Hedge's g)\nHeatwave - Control") + 
#   ylim(-5,5) + 
#   theme_pubr(base_size = 18)